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Abstract 

Recently, a new eigenvalue problem, called the transmission eigenvalue problem, has at¬ 
tracted many researchers. The problem arose in inverse scattering theory for inhomogeneous 
media and has important applications in a variety of inverse problems for target identification 
and nondestructive testing. The problem is numerically challenging because it is non-selfadjoint 
and nonlinear. In this paper, we propose a recursive integral method for computing transmis¬ 
sion eigenvalues from a finite element discretization of the continuous problem. The method, 
which overcomes some difficulties of existing methods, is based on eigenprojectors of compact 
operators. It is self-correcting, can separate nearby eigenvalues, and does not require an initial 
approximation based on some a priori spectral information. These features make the method 
well suited for the transmission eigenvalue problem whose spectrum is complicated. Numerical 
examples show that the method is effective and robust. 


1 Introduction 

The transmission eigenvalue problem (TE) [7, 4, 24, 5] has important applications in inverse scatter¬ 
ing theory for inhomogeneous media. The problem is non-selfadjoint and not covered by standard 
partial differential equation theory. Transmission eigenvalues have received significant attention in 
a variety of inverse problems for target identification and nondestructive testing since they provide 
information concerning physical properties of the target. 

Since 2010 significant effort has been focused on developing effective numerical methods for 
transmission eigenvalues [8, 25, 14, 28, 27, 1, 17, 6, 19]. The first numerical treatment appeared in 
[8], where three finite element methods were proposed. A mixed method (without a convergence 
proof) was developed in [14]. An and Shen [1] proposed an efficient spectral-element based numerical 
method for transmission eigenvalues of two-dimensional, radially-stratified media. The first method 
supported by a rigorous convergence analysis was introduced in [25]. In this article transmission 
eigenvalues are computed as roots of a nonlinear function whose values are eigenvalues of a related 
positive definite fourth order problem. This method has two drawbacks 1) only real transmission 
eigenvalues can be obtained, and 2) many fourth order eigenvalue problems need to be solved. In 
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[17] (see also [3]) surface integral and contour integral based methods are used to compute both real 
and complex transmission eigenvalues in the special case when the index of refraction is constant. 
Recently, Cakoni et.al. [6] reformulated the problem and proved convergence (based on Osborn’s 
compact operator theory [21]) of a mixed finite element method. Li et.al. [19] developed a finite 
element method based on writing the TE as a quadratic eigenvalue problem. Some non-traditional 
methods, including the linear sampling method in the inverse scattering theory [26] and the inside- 
out duality [18], were proposed to search for eigenvalues using scattering data. However, these 
methods seem to be computationally prohibitive since they rely on solving tremendous numbers 
direct problems. Other methods [10, 9, 13, 15] and the related source problem [11, 28] have been 
discussed in the literature. 

In general, developing effective finite element methods for transmission eigenvalues is challeng¬ 
ing because it is a quadratic, typically, degenerate, non-selfadjoint eigenvalue problem for a system 
of two second order partial differential equations and despite some qualitative estimates the spec¬ 
trum is largely unknown. In most cases, the continuous problem is degenerate with an infinite 
dimensional eigenspace associated with a zero eigenvalue. The system can be reduced to a single 
fourth order problem however conforming finite elements for such problems e.g. Argyris are ex¬ 
pensive. Straightforward finite element discretizations generate computationally challenging large 
sparse non-Hermitian matrix eigenvalue problems. Traditional methods such as shift and invert 
Arnoldi are handicapped by the lack of a priori eigenvalue estimates. To summarize, finite element 
discretizations of transmission eigenvalue problems generate large, sparse, typically highly degen¬ 
erate, non-Hermitian matrix eigenvalue problems with little a priori spectral information beyond 
the likelihood of a relatively high-dimensional nullspace. 

These characteristics suggest that most existing eigenvalue solver are unsuitable for transmission 
eigenvalues. Recently integral based methods [23, 22] related to the earlier work [12] and originally 
developed for electronic structure calculations become popular. These methods are based on eigen- 
projections [16] provided by contour integrals of the resolvent [2]. 

In this paper, we propose a recursive integral method (RIM) to compute transmission eigen¬ 
values from a continuous finite element discretization. Regions in the complex plane are searched 
for eigenvalues using approximate eigenprojections onto the eigenspace associated with the eigen¬ 
values within the region. The approximate eigenprojections are generated by approximating the 
resolvent contour integral around the boundary of the region by a quadrature on a random sample. 
The region is subdivided and subregions containing eigenvalues are recursively subdivided until the 
eigenvalues are localized to the desired tolerance. RIM is designed to approximate all eigenvalues 
within a specific region without resolving eigenvectors. This is well suited to the transmission eigen¬ 
value problem which typically seeks only the eigenvalues near but not at the origin. The degenerate 
non-hermitian nature of the matrix and the complicated unknown structure of the spectrum are 
not an issue. 

RIM is distinguished from other integral methods in literature by several features. First, the 
method works for Hermitian and non-Hermitian generalized eigenvalue problems such as those from 
the discretization of non-selfadjoint partial differential equations. Second, the recursive procedure 
automatically resolves eigenvalues near region boundaries and minimally separated eigenvalue pairs. 
Third, the method requires only linear solves with no need to explicitly form a matrix inverse. 

The paper is arranged as follows. Section 2 introduces the transmission eigenvalue problem, the 
finite element discretization, and the resulting large sparse non-Hermitian generalized matrix eigen- 
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value problem. Section 3 introduces the recursive integral method RIM to compute all eigenvalues 
within a region of the complex plane. Section 4 details various implementation details. Section 5 
contains results from a range of numerical examples. Section 6 contains discussion and future work. 

2 The transmission eigenvalue problem 

2.1 Formulation 

We introduce the transmission eigenvalue problem related to the Helmholtz equation. Let D CK ] 
be an open bounded domain with a Lipschitz boundary dD. Let k be the wave number of the 
incident wave u l = e lkx ' d and n{x) be the index of refraction. The direct scattering problem is to 
find the total field u(x) satisfying 


V • Vu + k 2 n(x)u = 0, 

in D , 

(la) 

A u + k 2 u = 0, 

in M 2 \ D, 

(lb) 

u(x) = e ikx - d + u s (x), 

on M 2 , 

(lc) 

- f dn s \ 

lim yft ( —- iku s ) = 0, 

r—> oo \ qt J 

where u s is the scattered field, x € M 2 ,r = \x\, d £ 

H := {x £ M 2 ; \x\ = 1}. 

(Id) 

The Sonmrerfeld 

radiation condition (Id) is assumed to hold uniformly with respect to x = x/\x\. 

The associated transmission eigenvalue problem is to find k such that there exist non-trivial 

solutions w and v satisfying 

V • Vw + k 2 n(x)w = 0, 

in D , 

(2a) 

Av + k 2 v = 0, 

in D , 

(2b) 

w — v = 0, 

on dD, 

(2c) 

dw dv 
~d^~d^^ ’ 

on dD, 

(2d) 


where u the unit outward normal to dD. The wave numbers L’s for which the transmission eigen¬ 
value problem has non-trivial solutions are called transmission eigenvalues. For existence results 
for transmission eigenvalues the reader is referred to the article and reference list of [5]. 

It is clear that k = 0 and w = v a harmonic function in D satisfies (2). So k = 0 is a non-trivial 
transmission eigenvalue with an infinite dimensional eigenspace. 

2.2 A continuous finite element method 

In the following, we describe a continuous finite element method for (2) [4, 13]. We use standard 
linear Lagrange finite element for discretization and define 

Sh = the space of continuous piecewise linear finite elements on D, 

s° h = s h n Hq(d) 

= the subspace of functions in Sh with vanishing DoF on dD, 

= the subspace of functions in Sh with vanishing DoF in D, 
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where DoF stands for degrees of freedom. 

Multiplying (2a) by a test function (p and integrating by parts gives 

(Vw, V4>) - k 2 (nw,(p) - = O’ (3) 

where (•,•) denotes the boundary integral on dD. Similarly, multiplying (2b) by a test function cp 
and integrating by parts gives 

(yv,V<P)-k 2 (v,<P)-(^,<p\ = 0. (4) 

Subtracting (4) from (3) and using the boundary condition (2d) gives 

(Vro — Vv, V<p) — k 2 ((nw — v), <p) = 0. (5) 

The Dirichlet boundary condition (2c) is explicitly enforced on the discretization by setting 

Wh = Wo,h + w B ,h where w 0;h E and w B , h € Sf, 

Vh = v 0>h + w B)h where v 0 ,h € 5°. 

Choosing the test function tph E S® for (3) gives the weak formulation for Wh as 

(VKh + WB,fi).V4) - k 2 {n(w 0th + w B)h ),£ h ) = 0, (6) 

for all £h E S°. Similarly, choosing the test function r//, € S'® gives the weak formulation for Vh as 

(y(vo,h + WB,h), V%) - k 2 ((v 0 , h + w B:h ),r] h ) = 0, (7) 

for all rjh € S Finally, choosing cp^ € S® in (5) gives 

(V(ic 0 ,/i + w B ' h ),V(p h ) - {'V(v 0 ,h + w B}h ),'V(ph) 

-k 2 {n(w 0th + WB,h) ~ (vo,h + WB,h), M = 0. (8) 

Let {£i,..., ^o} be the finite element basis for S° and 

be the basis for S),. Let N^. N°, and N® be the dimensions of Sh, S® and S&, respectively. Clearly 
{£jv 0 +i> • • • > £,N h } is a basis for S® and 

N h = N° + N*. 

Let S be the stiffness matrix given by (S)j j g = (V£j,V&), M n be the mass matrix given by 
(. M n )jj = (n£j,£e), and M be the mass matrix given by (M)jj = (£j,£g). Combining (6), (7), and 
(8), gives the generalized eigenvalue problem 

Ax = k 2 Bx, (9) 
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where matrices A and B are 


and 


( S N ° xN ° 


A = 


S N °h* N h 


S N ° xN h 

S N ° xN h 


\(S N h xN h) T (. -S N h xN h) T S N h xN h - S N h xN h 


( M, 


B = 


N ° xN h 


M N > N h 


, r N? x NP 

M n h h 
M n h xN h 


\ ( m5 xN *) t - ( M n ° xN h) t Mn* X N * - M N h X Ng 


A and B are clearly not symmetric and in general there are complex eigenvalues. Applications 
are typically interested in determining the structure of the spectrum (including complex conjugate 
pairs) near the origin. In practice, the primary focus is on computing a few of the non-trivial 
eigenvalues nearest the origin. Note for the transmission eigenvalue problem eigenvectors are of 
significantly less interest. 

Arnoldi iteration based adaptive search methods for real transmission eigenvalues were devel¬ 
oped in [14] and [20]. However, these methods are inefficient, may fail to converge, and are unable 
to compute all eigenvalues in general. The main goal of the current paper is to develop an effective 
tool to compute all the transmission eigenvalues (real and complex) of (9) in a region of the complex 
plane. 


3 A recursive contour integral method 

3.1 Continuous case 

We start by recalling some classical results in operator theory (see, e.g., [16]). Let T : X —> X be 
a compact operator on a complex Banach space X. The resolvent of T is defined as 

p(T) = {z € C : (z — T )~ x exists as a bounded operator on X}. (10) 

For any z € p(T), 

R Z (T) = (z- T)- 1 

is the resolvent of T and the spectrum of T is cr(T) = C \ p(T). 

Let T be a simple closed curve on the complex plane C lying in p(T) which contains m eigenvalues 
of T: Xi, i = 1,..., m. The spectral projection 

E[T] = hSr MT)dZ - 

is a projection onto the space of generalized eigenfunctions Ui,i = 1, ...,m associated with the 
eigenvalues A*,?' = 1,... ,m. If a function / has components in Ui,i = 1 then E(T)f is 

non-zero. If / has no components in u,, i = 1,..., m then E(T)f = 0. Thus E(T)f can be used to 
decide if a region contains eigenvalues of T or not. This is the basis of RIM. 

Our goal is to compute all the eigenvalues of T in a region S C C. RIM starts by defining 
T = dS, randomly choosing several functions fj,j = 1and approximating 

I j = E(T)f j , j = l,...,J, 


5 



by a suitable quadrature. Based on Ij we decide if there are eigenvalues inside S. If S contains 
eigenvalue(s), we partition S into subregions and recursively repeat this procedure for each subre¬ 
gion. The process terminates when each eigenvalue is isolated within a sufficiently small subregion. 

RIMO S,e,fjJ = l,...,J) 

Input: search region S, tolerance e, random functions fj,j = 1,... ,J 
Output: A, eigenvalue(s) of T in S 

1. Approximate (using a suitable quadrature) the integral 

E(T)f 3 = J R z {T)dzfj, j = r = 8S. 

2. Decide if S contains eigenvalue(s): 

— No. exit. 

— Yes. compute the size h(S) of S 

- if h(S) > e, partition S into subregions Si,i = 1,... N 

for i = 1 to N 

RIM (S i ,eJ j J = l,...,J) 

end 

- if h(S) < e, output the eigenvalue A and exit 

3.2 Discrete case 

We specialize RIM to potentially non-Hermitian generalized matrix eigenvalue problems. The 
finite element discretization of the transmission eigenvalue problem produces such a problem as do 
other similar discretizations of other PDEs. 

The matrix eigenvalue problem is 

Ax = A Bx (11) 

where A, B are n x n matrices, A is a scalar, and x is an n x 1 vector. The resolvent is 

R Z {A,B) = {zB- A)~ l B (12) 

for z in the resolvent set of the matrix pencil. The projection onto the generalized eigenspace 
corresponding to eigenvalues enclosed by a simple closed curve T is given by the Cauchy integral 

E(A, B) = — f {zB- A)~ 1 Bdz. (13) 

2m J T 

If the matrix pencil is non-defective then AX = BX A where A is a diagonal matrix of eigenvalues 
and X is an invertible matrix of generalized eigenvectors. This eigenvalue decomposition shows 

(zB - A)X = zBX - AX = zBX - BXA = BX(zI - A), 
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and gives 


X(zl - A) _i Y _i = (zB - A)~ l B 

for complex 2 : not equal to any of the eigenvalues. Integrating the resolvent around a closed contour 
r in C gives 


1 

2iri 


R Z (A, B)dz = X 


1 

2tti 


(zl - A )- l dzX- ^ = XK V X 


-1 


—1 


where Ap is A with eigenvalues inside T set to 1 and those outside T set to 0. 

The projection of a vector y onto the generalized eigenspace for eigenvalues inside T is 

1 


Py := XhrX-'y = 


2ni 


R Z (A, B)ydz. 


(14) 


If there no eigenvalues are inside T, then P = 0 and Py = 0 for all y £ C n . 
We select a quadrature rule to approximate the contour integral 


N 

<?=! 


where ui q and z q are the quadrature weights and points, respectively. Although an explicit compu¬ 
tation of R z is not possible one can approximate the projection of y by 


N 

<7=1 


where r q are the solutions of the linear systems 

1 


(■ z q B - A)r q = —uj q By, q = 1,. 


, N. 


(15) 


(16) 


For robustness, we use a set of vectors y 3 , j = 1,..., J assembled as the columns of an n x J 
matrix Y. The RIM for generalized eigenvalue problems is as follows. 

M-RIM (A,B,S,e,Y) 

Input: matrices A and B , search region S, tolerance e, random vectors Y 
Output: generalized eigenvalue A 

1. Compute Py 3 , j = 1,..., J using (15) on dS. 

2. Decide if S contains eigenvalue (s): 

— No. exit. 

— Yes. compute the size h(S) of S 

- if h(S) > e, partition S into subregions Si,i = 1,... I 

for i = 1 to I 

M-RIM (A,B,Si,e,Y) 

end 

- if h(S) < e, output the eigenvalue A and exit 
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4 Implementation 


We assume the search region S' is a polygon in the complex plane C for simplicity and divide S 
into subregions of simple geometry, such as triangles and rectangles. Rectangles are used in the 
implementation. 

There are several keys in the implementation of RIM: we need a suitable quadrature rule for 
the contour integral; we need a mechanism to solve (16); and we need an effective rule to decide if 
a subregion contains eigenvalues. 

We use Gaussian quadrature on each rectangle edge. It does not appear necessary to use many 
points and we use the two point rule. 

In contrast with the quadrature, an accurate linear solver seems necessary and we use the 
Matlab “\” command. 

Next we discuss the rule to decide if S might contain eigenvalues and needs to be subdivided. 
We refer to a subregion that potentially contains at least one eigenvalue as admissible. Any vector y 
is represented in the eigenbasis (columns of A) as y = a i x i- Assume there are M eigenvalues 
inside T and reorder the eigenvalues and eigenvectors with these M eigenvectors as xi, X 2 ,..., xy 
then 

M 

Py = — R z (A,B)ydz = X\ r X 1 y = VV x i . 

2m J r tT 

So it is reasonable to use ||Py|| to decide if a region contains eigenvalues. There are two primary 
concerns for the robustness of the algorithm. We might miss eigenvalues if ||Py|| is small when 
there is an eigenvalue within T. We might continue to subdivide a region if | Py \ | is large when 
there is no eigenvalue within T. In the first case ||Py|| could be small when there is an eigenvalue 
because of quadrature/rounding errors and/or simply because the random components a* are small. 
Our solution is to project Py again with an amplifier K and look at ||P(IFPy)||. In fact, one can 
simply choose K = l/||Py||. In the second case ||Py|| could be large when there is no eigenvalue 
inside T if there are eigenvalues right outside T and the quadrature rule or the linear solver are not 
sufficiently accurate. Fortunately, RIM has an interesting self-correction property that fixes such 
errors on subsequent iterations. 

In our implementation, we use the following rules to decide an admissible region: 

1. We use several random vectors y j,j = 1,..., J; 

2. We use |P(RPyj) || where K is an amplifier. 

Rule 1. and Rule 2. guarantee that even if the component of y in X is small, the algorithm can 
detect it effectively since P(KPy) should be of the same size of KP{y). If there is no eigenvalue 
inside T, Py can still be large due to reasons we mentioned above. However, another projection of 
P(KPy) should significantly reduce ||/iPy||. 

The indicator function xs is the ratio 

, JjP(gPy)j! 

xs ■ \\Py\\ ■ 
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\\P(KPy)\\ 

If there are eigenvalues inside T, then -——— = 0(K). On the contrary, if there is no eigenvalue 


i-ide r, = «*>. 


l|J=y|l 


Here are some details in the actual implementation. 

1. The search region S' is a rectangle. 

2. We use 3 random vectors y j,j = 1, 2,3. 

3. The amplifier is set as K = 10. 

4. We use 2 point Gauss quadrature rule on each edge of S. 

5. We use Matlab ”\” to solve the linear systems. 

6. We take the indicator function as 

\\P(KPyj))\ 


XS = max 
3 = 1,2,3 


\\p yj 


7. We use K /10 as the criterion, i.e., if xs > Al/10, S is admissible. 


5 Numerical Examples 

In this section, we assume that the initial search region S is a rectangle. We present examples to 
show the performance of RIM. 

5.1 Transmission Eigenvalues 

We test RIM on the generalized matrix eigenvalue problem for transmission eigenvalues using 
continuous finite element method described in Section 2. Since the original partial differential 
problem is non-selfadjoint, the generalized matrix eigenvalue problem is non-Hermitian. In practice, 
we only need a few eigenvalues of smallest norm. However, we do not have an a prior knowledge of 
the locations of the eigenvalues. 

Example 1: We consider a disc D with radius 1/2 and index of refraction n(x) = 16 where 
the exact transmission eigenvalues [8] are the roots of 

Ji(k/2)Jo(2k) = 4Jo(fc/2) J±(2k), m = 0, 

and 

Jm-i(k/2)J m (2k) = 4J m (k/2)J m _i(2k), m > 1, 
where J m ’s are Bessel functions. 

A regular mesh with with h ~ 0.05 is used to generate the 1018 x 1018 matrices A and B and 
we consider the preliminary search region S = [1,10] x [—1,1], Since the mesh is relatively coarse 
we take e = l.Oe — 3 and use 3 random vectors. RIM computes 3 eigenvalues 

Ai = 3.994, A 2 = 6.935, A 3 = 6.939 
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which are good approximations of the exact eigenvalues given in [4] 

Ai = 3.952, A 2 = 6.827, A 3 = 6.827. 

Note that the values we compute are k 2, s and the actual values in [4] are k’s. 

As a second test we choose S = [22,25] x [—8,8] and find 4 eigenvalues in this region 

Ax = 24.158 + 5.690i, A 2 = 24.158 - 5.690i, A 3 = 25.749, A 4 = 25.692 

which approximate the exact eigenvalues 

Ai i2 = 23.686 ± 5.667i, A 3)4 = 24.465. 

Note that RIM computes the generalized eigenvalues to the anticipated accuracy e the discrepancy 
is mainly due to the fact finite element methods approximate smaller eigenvalues better than larger 
eigenvalues. 

The search regions for the transmission eigenvalue tests are shown in Fig 1. The algorithm 
refines near the eigenvalues until the tolerance is met. The right image in Fig. 1 shows only three 
refined regions because two eigenvalues are very close. 




Figure 1: The regions explored by RIM for the disc with radius 1/2, n{x) = 16, and e = l.Oe — 3. 
Left: the search region is given by S = [1,10] x [—1,1]. Right: the search region is given by 
S = [22,25] x [-8,8]. 

Example 2: Let D be the unit square and n(x) = 16 (AAS: I think this is correct) with 
h ~ 0.05. The matrices A and B are 2075 x 2075. The exact transmission eigenvalues are not 
available. The first search region is given by S' = [3,8] x [—1,1]. RIM computes the following 
eigenvalues 

Ai = 3.561, A 2 = 6.049, A 3 = 6.051. 

They are consistent with the values given in Table 3 of [4]: 

Ai = 3.479, A 2 = 5.883, A 3 = 5.891. 
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The second search region is given by S = [20,25] x [—8,8]. The eigenvalues we obtain are 

Ai = 20.574 + 5.128*, A 2 = 20.574 - 5.128*, A 3 = 21.595, A 4 = 23.412. 

We plot the search regions in Fig. 2. The left picture is for S = [3, 7] x [—1,1], The right picture 
is for S = [20,25] x [-6,8], 





F=~ 


a 


22.5 

real axis 


Figure 2: The regions explored by RIM for the unit square with n{x) = 16 and e = l.Oe — 3. Left: 
the search region is given by S = [3, 7] x [—1,1], Right: the search region is given by [20, 25] x [—6, 8]. 


5.2 Eigenvalues on T := dS 

It is very unlikely that T := dS is not contained in the resolvent set. However, we want to explore 
what will happen if eigenvalues he on on T. The first example shows that this does not generate 
difficulty for RIM. 

Example 3: We first consider a simple example given below (Example 5 of [23]): 


/ 99 
100 
0 


A = 



100 


0 

V o 


0 

0 

0 


0 \ 
0 


1 

100 


0 


r 

100/ 


B = diag(0,..., 0,1,..., 1), 


where B has 20 ones on its diagonal. The following are some exact eigenvalues 


Ai =0, A 2 = 0.01, A 3 = 0.02, A 4 = 0.03. 

We set the initial search region to be S = [0,1/30] x [0,1/100] and e = l.Oe — 9 and note that 
all the eigenvalues are on T := dS. 
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real axis 


Figure 3: Eigenvalues on F = dS. All the four eigenvalues are on T. 


The eigenvalues computed by RIM are given below (see also Fig. 3). They are accurate up to 
the required precision. From Fig. 3, we can see that RIM keeps refining around the eigenvalues. 

Ai = (4.967053731282552 + 4.967053731282552*)10 -10 , 

A 2 = 0.009999999900659 + 0.000000000496705i, 

A 3 = 0.020000000298023 + 0.000000000496705i, 

A 4 = 0.029999999701977 + 0.000000000496705h 


5.3 Self-correction Property 

When a quadrature point z q in the collection of linear systems (16) is close to an eigenvalue A, the 
linear system will be ill-conditioned. In particular, when A is just outside T the indicator function 
Xs could be large because either the linear solve or quadrature rule are not sufficiently accurate. 
RIM will take such regions as admissible and refine. But fortunately, after a few subdivisions, 
RIM appears to discard the sub regions. We demonstrate this interesting self-correction property 
using two example. 

Example 4: We use matrices A and B from Example 2 and focus on the eigenvalue located 
at 3.9945. We choose the initial search region S = [4.0,4.2] x [0, 0.2] and note that there is no 
eigenvalue in S. With the same standard two-point Gauss quadrature rule on each edge of S RIM 
computes 

XS = 4.383, (17) 
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indicating that there may be eigenvalues in S and RIM procedes to recursively explore S by 
dividing S into the four rectangles 

S\ = [4.0,4.1] x [-0.1,0], Si = [4.0,4.1] x [0,0.1], 

S 3 1 = [4.1,4.2] x [-0.1,0], Si = [4.1,4.2] x [0,0.1] 

with indicator values 

Xs\ = L589 > X.S’1 = 1-589, Xsi = 0.002, xs\ = °- 002 - 

RIM discards S 3 and S\ and retains Sj and as admissible regions. 

We show the result for region S{: Si is similar. The four rectangles from dividing Sj are 

Sf = [4.0,4.05] x [-0.05,0], S’| = [4.0,4.05] x [-0.1,-0.05], 

S‘i = [4.05,4.1] x [-0.05,0], Sj = [4.0,4.05] x [-0.1,-0.05], 

with indicator values 

X S f = 0.997, xsi = °- 002 > XSg = °- 002 i Xsi = 2 - 159e “ ° 4 - 

and RIM discards all the regions. Let us see one more level. Suppose Xsf an4 subdivided 

into 


Sf = [4.0,4.025] x [-0.025,0], S% = [4.0,4.025] x [-0.05,-0.025], 

Sf = [4.025,4.05] x [-0.025,0], Sf = [4.025,4.05] x [-0.05,-0.025] 

with indicator values 

Xsi = °- 395 > Xsi = °- 002, Xs l = °- 001, Xs l = 1M5e ~ ° 4 - 
Hence RIM eventually discards S. 

Example 5: The same experiment is conducted for a search region around the complex eigen¬ 
value A = 24.1586 + 5.690?' with initial search region S = [24.16,24.96] x [5.30,6.10] which although 
close to the eigenvalue does not contain any eigenvalues. Indicator values are in Table. 1 and we 
can note that RIM does eventually conclude that there are no eigenvalues in the region. 


5.4 Close Eigenvalues 

RIM is able to separate nearby eigenvalues provided the tolerance is less than the eigenvalue 
separation. 

Example 6: This example comes from a finite element discretization of the Neumann eigenvalue 
problem: 


—A u = Ait, 
dii 


in H, 

(18a) 

on ()I). 

(18b) 
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Table 1: The indicators function \S on different search regions. 


[24.16.24.56] x [5.30,5.70] 11.825 

[24.56,24.96] x [5.30,5.70] 5.418e-ll 

[24.16.24.36] x [5.30,5.50] 9.216e-ll 

[24.36.24.56] x [5.30,5.50] 8.712e-14 

[24.16.24.26] x [5.50,5.60] 1.742e-ll 

[24.26.24.36] x [5.50,5.60] 1.476e-13 

[24.16.24.21] x [5.60,5.65] 6.558e-10 

[24.21.24.26] x [5.60,5.65] 1.378e-13 

[24.16,24.185] x [5.65,5.675] 1.159e-8 

[24.185.24.21] x [5.65,5.675] 4.000e-13 

[24.16,24.185] x [5.65,5.675] 5.574e-06 

[24.185,24.21] x [5.65,5.675] 4.304e-12 


[24.16.24.56] x [5.70,6.10] 0.195 

[24.56,24.96] x [5.70,6.10] 4.119e-ll 

[24.16.24.36] x [5.50,5.70] 3.682 

[24.36.24.56] x [5.50,5.70] 5.870e-ll 

[24.16.24.26] x [5.60,5.70] 7.806 

[24.26.24.36] x [5.60,5.70] 6.755e-ll 

[24.16.24.21] x [5.65,5.70] 2.799 

[24.21.24.26] x [5.65,5.70] 8.229e-ll 

[24.16,24.185] x [5.675,5.70] 1.556 

[24.185.24.21] x [5.675,5.70] 8.648e-ll 

[24.16,24.1725] x [5.6875,5.70] 0.095 

[24.185,24.21] x [5.675,5.70] 2.628e-ll 


where D is the unit square which has an eigenvalue 7r 2 of multiplicity 2. We use linear Lagrange 
elements on a triangular mesh with h ~ 0.025 to discretize and obtain a generalized eigenvalue 
problem 

Ax = ARx, (19) 

where the stiffness matrix A and mass matrix B are 2075 x 2075. The discretization has broken 
the symmetry and (19) the eigenvalue of multiplicity 2 has been approximated by a very close pair 
of eigenvalues of 

Ax = 9.872899741642826 and A 2 = 9.872783160389966. 

With e = l.Oe — 3 RIM fails to separate the eigenvalues and we obtain only one eigenvalue 

A! = 9.872680664062500. 

However, with e = l.Oe — 9 RIM separates the eigenvalues and we obtain 

Ax = 9.872899741516449 and A 2 = 9.872783160419203. 

The search regions explored by RIM with different tolerances are shown in Fig. 4. 

Example 7: As a final example we compute the eigenvalues of the 40 x 40 Wilkinson matrix 
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x 10~ 4 




Figure 4: The regions explored by RIM. The search region is given by S = [1,10] x [—1,1]. Left: 
e = l.Oe — 3. Right: e = l.Oe — 9. 

which is known to have very close eigenvalues. With e = l.Oe — 14 and the search region S = 
[-2,10] x [-2, 10]. RIM accurately distinguishes the close eigenvalues with giving the results 
shown in Table 2 and Fig. 5. 


Table 2: The computed Wilkinson eigenvalues by RIM. 


i 

A* 

i 

Ai 

1 

-1.125441522046458 

11 

5.000236265619321 

2 

0.253805817279499 

12 

5.999991841327017 

3 

0.947534367500339 

13 

6.000008352188331 

4 

1.789321352320258 

14 

6.999999794929806 

5 

2.130209219467361 

15 

7.000000207904748 

6 

2.961058880959172 

16 

7.999999996191775 

7 

3.043099288071971 

17 

8.000000003841876 

8 

3.996047997334983 

18 

8.999999999945373 

9 

4.004353817323874 

19 

9.000000000054399 

10 

4.999774319815003 

20 

9.999999999999261 


6 Discussion and future work 

This paper proposes a robust recursive integral method RIM to compute transmission eigenvalues. 
The method effectively locates all eigenvalues in a region when neither the location or number 
eigenvalues is known. The key difference between RIM and other counter integral based methods 
in the literature is that RIM essentially only tests if a region contains eigenvalues or not. As a 
result accuracy requirements on quadrature, linear solves, and the number of test vectors may be 
significantly reduced. 
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Figure 5: The regions explored by RIM for the Wilkinson matrix with e = l.Oe — 14. 


RIM is a non-classical eigenvalue solver which is well suited to problems that only require 
eigenvalues. In particular, the method snot only works for matrix eigenvalue problems resulting 
from suitable numerical approximations, e.g., finite element methods, of PDE-based eigenvalue 
problem, but also those eigenvalue problems which can not be easily casted as a matrix eigenvalue 
problem, e.g., see [3, 17]. 

The goal of this paper is to introduce the idea of RIM and demonstrate its potential to compute 
eigenvalues. A paper like this raises more questions than it answers. How inaccurate can the 
quadrature be and still locate eigenvalues? How inaccurate can the the linear solver can and still 
locate eigenvalues. The current implementation uses a combination of inaccurate quadrature and 
accurate solver: two point Gaussian quadrature on the edges of rectangles and the Matlab “\” 
operator. These two separate issues can be combined into one question: how accurate does the 
overall procedure have to be to accurately distinguish admissible regions. These crucial complexity 
issues are not addressed in this current paper. 

The example problems are small. We plan to extend RIM for large (sparse) eigenvalue problems 
which will require replacing “\” with an iterative solver. Parallel extension is another interesting 
project since the algorithm is essentially embarrassingly parallel. In particular, a GPU implement 
of RIM is under consideration. 
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